#include "mt3d.h"

void mt3d::apply_boundary_condition()
{
    std::clog << "Applying boundary condition ... ";

    // 采用划0置1的方法处理第一类边界条件 （稀疏矩阵处理起来很快）

    for (int i = 0; i < edge_num_; ++i)
    {
        if (edge_on_boundary_[i])
        {
            B_TE_.coeffRef(i) = _zero; // 二次场方法计算中假设边界上的二次场已经衰减为0
            B_TM_.coeffRef(i) = _zero;
        }
    }

    for (int i = 0; i < edge_num_; ++i) // 外层是行循环
    {
        if (!edge_on_boundary_[i])
        {
            for (Eigen::SparseMatrix<complex_d, Eigen::RowMajor>::InnerIterator it(A_PDE_, i); it; ++it) // 列循环
            {
                if (edge_on_boundary_[it.col()])
                {
                    B_TE_[i] -= it.value() * B_TE_[it.col()];
                    B_TM_[i] -= it.value() * B_TM_[it.col()];
                }
            }
        }
    }

    for (int i = 0; i < edge_num_; ++i) // 处理列
    {
        for (Eigen::SparseMatrix<complex_d, Eigen::RowMajor>::InnerIterator it(A_PDE_, i); it; ++it)
        {
            if (edge_on_boundary_[it.col()])
            {
                it.valueRef() = _zero;
            }
        }
    }

    for (int i = 0; i < edge_num_; ++i) // 处理行
    {
        if (edge_on_boundary_[i])
        {
            for (Eigen::SparseMatrix<complex_d,  Eigen::RowMajor>::InnerIterator it(A_PDE_, i); it; ++it)
            {
                it.valueRef() = _zero;
                if (it.row() == it.col()) it.valueRef() = _one; // 置1
            }
        }
    }

    A_PDE_.makeCompressed();

    std::clog << "done\n";
    return;
}
